/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 * Marcus Cuda
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
namespace dnAnalytics.LinearAlgebra.Decomposition
{
    internal class DenseLU : AbstractLU
    {
        public DenseLU(Matrix matrix) : base(matrix)
        {
        }

        protected override void DoCompute(Matrix data, int[] pivots)
        {

            DenseMatrix matrix = (DenseMatrix) data;

            int order = matrix.Rows;
            for (int i = 0; i < order; i++)
            {
                pivots[i] = i;
            }

            double[] LUcolj = new double[order];

            // Outer loop.
            for (int j = 0; j < order; j++)
            {
                int indexj = j*order;
                int indexjj = indexj + j;
                // Make a copy of the j-th column to localize references.
                for (int i = 0; i < order; i++)
                {
                    LUcolj[i] = matrix.Data[indexj + i];
                }

                // Apply previous transformations.
                for (int i = 0; i < order; i++)
                {
                    // Most of the time is spent in the following dot product.
                    int kmax = System.Math.Min(i, j);
                    double s = 0.0;
                    for (int k = 0; k < kmax; k++)
                    {
                        s += matrix.Data[k*order + i]*LUcolj[k];
                    }

                    matrix.Data[indexj + i] = LUcolj[i] -= s;
                }

                // Find pivot and exchange if necessary.
                int p = j;
                for (int i = j + 1; i < order; i++)
                {
                    if (System.Math.Abs(LUcolj[i]) > System.Math.Abs(LUcolj[p]))
                    {
                        p = i;
                    }
                }
                if (p != j)
                {
                    for (int k = 0; k < order; k++)
                    {
                        int indexk = k*order;
                        int indexkp = indexk + p;
                        int indexkj = indexk + j;
                        double temp = matrix.Data[indexkp];
                        matrix.Data[indexkp] = matrix.Data[indexkj];
                        matrix.Data[indexkj] = temp;
                    }
                    pivots[j] = p;
                }

                // Compute multipliers.

                if (j < order & matrix.Data[indexjj] != 0.0)
                {
                    for (int i = j + 1; i < order; i++)
                    {
                        matrix.Data[indexj + i] /= matrix.Data[indexjj];
                    }
                }
            }
        }

        protected override Matrix ComputeInverse(Matrix matrix, int[] pivots)
        {
            int order = matrix.Rows;
            DenseMatrix inverse = new DenseMatrix(order, order);
            for (int i = 0; i < order; i++)
            {
                inverse.Data[i + order*i] = 1.0;
            }

            Solve(mMatrix, mPivots, inverse);
            return inverse;
        }


        private static void Pivot(int m, int n, double[] B, int[] pivots)
        {
            for (int i = 0; i < pivots.Length; i++)
            {
                if (pivots[i] != i)
                {
                    int p = pivots[i];
                    for (int j = 0; j < n; j++)
                    {
                        int indexk = j*m;
                        int indexkp = indexk + p;
                        int indexkj = indexk + i;
                        double temp = B[indexkp];
                        B[indexkp] = B[indexkj];
                        B[indexkj] = temp;
                    }
                }
            }
        }

        private static void Solve(int order, int columns, double[] factor, int[] pivots, double[] data)
        {
            Pivot(order, columns, data, pivots);

            // Solve L*Y = B(piv,:)
            for (int k = 0; k < order; k++)
            {
                int korder = k*order;
                for (int i = k + 1; i < order; i++)
                {
                    for (int j = 0; j < columns; j++)
                    {
                        int index = j*order;
                        data[i + index] -= data[k + index]*factor[i + korder];
                    }
                }
            }
            // Solve U*X = Y;
            for (int k = order - 1; k >= 0; k--)
            {
                int korder = k + k*order;
                for (int j = 0; j < columns; j++)
                {
                    data[k + j*order] /= factor[korder];
                }
                korder = k*order;
                for (int i = 0; i < k; i++)
                {
                    for (int j = 0; j < columns; j++)
                    {
                        int index = j*order;
                        data[i + index] -= data[k + index]*factor[i + korder];
                    }
                }
            }
        }

        protected override void Solve(Matrix factor, int[] pivots, Matrix result)
        {
            Solve(factor.Rows, result.Columns, ((DenseMatrix) factor).Data, pivots, ((DenseMatrix) result).Data);
        }

        protected override void Solve(Matrix factor, int[] pivots, Vector result)
        {
            Solve(factor.Rows, 1, ((DenseMatrix) factor).Data, pivots, ((DenseVector) result).Data);
        }
    }
}